% Specify the folder containing the images
folderPath = ''; % Change this to your folder path

% Get all PNG image files
imageFiles = dir(fullfile(folderPath, '*.jpg'));

% Check if any images exist
numImages = length(imageFiles);
if numImages == 0
    disp('No images found in the specified folder.');
    return;
end

% Create a subfolder for saving the results
resultsFolder = fullfile(folderPath, 'Filtered_Clusters');
if ~exist(resultsFolder, 'dir')
    mkdir(resultsFolder);
end

% Set the minimum cluster area threshold (adjust as needed)
minClusterSize = ; % Change this value to filter smaller or larger clusters

% Initialize Excel file path
excelFilePath = fullfile(resultsFolder, 'Filtered_Cluster_Properties.xlsx');

% Loop through all images in the folder
for i = 1:numImages
    try
        % Read the already binarized and filled image
        imageName = imageFiles(i).name;
        imagePath = fullfile(folderPath, imageName);
        bw_clean = imread(imagePath);

        % Ensure image is grayscale (some binarized images might still be RGB)
        if size(bw_clean, 3) == 3
            bw_clean = rgb2gray(bw_clean);
        end

        % Convert to logical format if it's not already
        bw_clean = bw_clean > 0;

        % Remove small noise from the binary image (preliminary filtering)
        bw_clean = bwareaopen(bw_clean, minClusterSize);

        % Invert image so clusters are labeled correctly (black = foreground, white = background)
        bw_inverted = ~bw_clean;  

        % Label connected components
        labeledImage = bwlabel(bw_inverted);

        % Get properties of the clusters (Area, Centroid, Perimeter)
        stats = regionprops(labeledImage, 'Area', 'Centroid', 'Perimeter', 'PixelIdxList');

        % **Filter clusters based on area threshold**
        validClusters = [];
        for j = 1:length(stats)
            if stats(j).Area >= minClusterSize
                validClusters = [validClusters; j]; % Store valid cluster indices
            else
                labeledImage(stats(j).PixelIdxList) = 0; % Remove small clusters
            end
        end

        % Relabel the remaining clusters
        labeledImage = bwlabel(labeledImage > 0);

        % **Generate a unique colormap for the remaining clusters**
        numClusters = max(labeledImage(:));
        uniqueColors = distinguishable_colors(numClusters, [1 1 1]); % Avoid white background

        % Convert labeled image to RGB using the unique color map
        filteredImageColor = label2rgb(labeledImage, uniqueColors, 'w', 'shuffle');

        % Get the centroid coordinates of valid clusters
        centroids = cat(1, stats(validClusters).Centroid); % Extract centroids
        perimeters = [stats(validClusters).Perimeter]'; % Extract perimeters
        clusterIDs = (1:length(validClusters)); % Cluster ID numbers

        % **Overlay cluster numbers using text()**
        figure;
        imshow(filteredImageColor);
        hold on;
        for k = 1:length(clusterIDs)
            text(centroids(k,1), centroids(k,2), num2str(clusterIDs(k)), ...
                'Color', 'black', 'FontSize', 14, 'FontWeight', 'bold', ...
                'HorizontalAlignment', 'center');
        end
        hold off;

        % Generate filenames for output
        [~, fileName, ext] = fileparts(imageFiles(i).name);
        outputImagePath = fullfile(resultsFolder, [fileName '_filtered_clusters.png']);

        % Save the annotated image
        saveas(gcf, outputImagePath);
        close; % Close figure to avoid excessive open windows

        % Confirm file save
        if exist(outputImagePath, 'file')
            disp(['Successfully saved: ', outputImagePath]);
        else
            warning(['Failed to save: ', outputImagePath]);
        end

        % Save the filtered cluster properties to an Excel sheet
        if ~isempty(validClusters)  % Ensure there are clusters to save
            % Extract only the valid cluster properties
            filteredStats = stats(validClusters);

            % Save data if any valid clusters remain
            if ~isempty(filteredStats)
                clusterData = table((1:length(filteredStats))', [filteredStats.Area]', ...
                                    centroids(:,1), centroids(:,2), perimeters, ...
                                    'VariableNames', {'ClusterID', 'Area', 'Centroid_X', 'Centroid_Y', 'Perimeter'});

                % Excel sheet names must be <=31 characters
                sheetName = fileName(1:min(31, length(fileName)));

                % Save without overwriting the whole file
                writetable(clusterData, excelFilePath, 'Sheet', sheetName, 'UseExcel', false);
            end
        end

        % Display progress message
        disp(['Processed and saved: ', outputImagePath]);

    catch ME
        % Display an error message if something goes wrong
        warning(['Error processing ', imageFiles(i).name, ': ', ME.message]);
    end
end

disp('All images have been processed. Filtered clusters with IDs are saved in the "Filtered_Clusters" folder.');

% Function to generate distinguishable colors
function colors = distinguishable_colors(n_colors, bg)
    if nargin < 2
        bg = [1 1 1]; % Default white background
    end
    colors = lines(n_colors); % Default MATLAB colormap (modify for larger color diversity)
    if n_colors > size(colors, 1)
        colors = rand(n_colors, 3); % Generate extra colors if needed
    end
end
